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Techniques for simulating molecules whose conformations satisfy constraints are presented. A method for select- 
ing appropriate moves in Monte Carlo simulations is given. The resulting moves not only obey the constraints 
but also maintain detailed balance so that correct equilibrium averages are computed. In addition, techniques for 
optimizing the evaluation of implicit solvent terms are given. 
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I. INTRODUCTION 

When attempting to compute thermodynamic quantities 
with a molecular simulation, we are frequently confronted 
with the problem of sampling in a high-dimensional config- 
uration space. The dimensionality of this space is given by 
the number of degrees of freedom for the molecular system. 
Techniques which lower the number of degrees of freedom 
will increase the efficiency of the thermodynamic sampling — 
provided, of course, that these techniques are physically justi- 
fied. Thus, an implicit solvent model may be used to eliminate 
the degrees of freedom associated with the solvent molecules; 
the standard chemical force fields replace the electron charges 
with atom-centered partial charges thereby removing the elec- 
trons' degrees of freedom. Further reductions in dimensional- 
ity are possible by imposing constraints on the relative posi- 
tions of the atoms in a molecule. Thus we might specify that 
the bond lengths and bond angles in a molecule are fixed and 
only the torsion angles are allowed to vary. It is such a sce- 
nario that we examine in this paper. We address two aspects of 
this problem: how to move a molecule subject to constraints 
in order to allow equilibrium averages to be computed using 
the canonical-ensemble Monte Carlo method 1 1 ] and how to 
evaluate the energy efficiently. 

The imposition of constraints in molecular modeling has 
been extensively studied |2, §3.3.2, §15.1]. Let us start by elu- 
cidating the difference in the treatment of hard constraints in 
molecular dynamics and Monte Carlo simulations. We treat 
hard constraints by taking the limit where the "spring con- 
stant" for the hard degrees of freedom is infinite. Molecular 
dynamics simulations then consider the evolution of the re- 
sulting system over a finite time. On the other hand, if we 
wish to determine the equilibrium properties of a system using 
the Monte Carlo method, then we need to consider averaging 
over sufficiently long times to allow equipartition of energy 
among all the degrees of freedom of a system. This is, of 



course, an example of nonuniform limits. We are interested in 
taking both t s „ 



oo and T cqu — > oo, where T s ; m is the rep- 
oqu is the equipartition time 



resentative simulation time and r, 
(which is proportional to the "stiffness" of the constraints). In 
constrained molecular dynamics, we take the limit r cqu — ► oo 
first, which prevents equipartition from occurring; whereas in 
equilibrium statistical mechanics, we take T s ; m — > oo first and 
this allows energy equipartition. If we are attempting to com- 
pute an equilibrium quantity, such as the free energy of bind- 
ing, it is essential to allow energy equipartition. Understand- 
ing this distinction explains the apparently contradictory re- 
sults for constrained and unconstrained averages for a flexible 
ti-imer|2, §15.1]. 

One way of understanding the constrained equilibrium sys- 
tem is to consider how the equilibrium distribution varies as 
the constraint is imposed. In the limit, the distribution col- 
lapses to a lower-dimensional sub-manifold of configuration 
space. However, this sub-manifold has a "thickness" that de- 
pends on the details of the constraint term and, consequently, 
Monte Carlo moves for the constrained system need to reflect 
this thickness in order to sample the distribution correctly. As 
a consequence, we will need to specify the functional form of 
the constraint energy and the constraint is no longer a purely 
geometrical object. At first glance, this would appear to com- 
plicate further the already complex algebra of constrained mo- 
tions 01 • However, we will propose an algorithm for making 
moves which is simple to implement and which automatically 
ensures that the correct equilibrium averages are computed. 

The second half of the paper considers a mundane — but 
nevertheless important — problem, namely how to evaluate the 
energy of a molecule made up of rigid subcomponents. We 
propose a consistent framework for avoiding the computation 
of constant terms and for imposing energy cutoffs. We extend 
this to the computation of the generalized Born solvation term 
and we describe a simple method for computing the solvent 
accessible surface area which has a bounded error. 
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II. GENERALIZED MONTE CARLO MOVES 

We begin by assembling some techniques for combining 
Monte Carlo moves. We define an "E move" as an ergodic 
move which preserves exp(— j3E) as the invariant distribu- 
tion, where (3 = l/(kT) and k is the Boltzmann constant. 
(Here, "ergodic" implies that the move allows all relevant 
portions to configuration space to be explored.) Note that 
a zero move has a uniform invariant distribution. A typi- 
cal zero move samples a new configuration from a distribu- 
tion which satisfies the symmetry requirement p(T';T) = 
p(T; r'), where p(T'; T) is the probability density of picking 
a new configuration of V given a starting configuration of T. 
Clearly a sequence of n E moves is itself an E move. From 
the central limit theorem, a sequence of n zero moves is equiv- 
alent, in the limit of large n, to selecting the new configuration 
from a multi-dimensional Gaussian. 

Instead of carrying out the n moves with a given energy 
E(T), we can consider the case where the energy is given by 
E\(T) which depends continuously on the parameter A. A 
sequence of n E\ moves where A is varied adiabatically in 
such a way that its initial and final values are Ao is an E\ 
move. This follows because adiabatically varied systems are 
always in equilibrium with the instantaneous value of A 0. 
§11]. Each E\ move is carried out at a fixed A and A is varied 
between the moves. In order to satisfy the adiabatic condition, 
we will need to take n large. 

A move from T to V may be subjected to "Boltzmann ac- 
ceptance with energy E". This involves accepting the move 
(V is the new state) with probability M(x) and otherwise 
rejecting the move (T is the new state). Here M(x) is a 
function satisfying < M(x) < 1 and M{x)/M(—x) — 
exp(-x) with x = 0(E(T') - E(T)). Usually we take 
M(x) = min(l, cxp(— x)); however other choices, e.g^, the 
Fermi function, M(x) = 1/(1 + exp(x)), are possible |5]. 

Consider an E\ move from T to V followed by an Boltz- 
mann acceptance using E2- This compound move is an 
(Ei + E2) move. The proof follows as a special case of 
the "multiple time-step" (MTS) method \&l §11] or, alterna- 
tively, as a special case of early rejection |2|, §14.3.2]. If the 
Ei move was already a "rejected" move, i.e., V = T, then 
the Boltzmann test involving E2 automatically "succeeds" 
(M (0) = 1). Thus E2 does not need to be evaluated in this 
case. 

These results allow us to generalize the MTS method by 
splitting the energy into to terms (instead of just two), 

m 

£7(T) = £l5,(r). 
z=i 

The method is defined recursively as follows: a level-0 move 
is defined to be a zero move; a level-Z move, with I > 0, is 
defined to be n;_i level- (I — 1) moves the result of which is 
subjected to Boltzmann acceptance using E\. By induction, 
we see that a level-Z move is an £\ move, where 

1 

£(r) = £^,(r). 
i'=i 



It follows that a level-m move is an £ m move, i.e., an E move. 
Typically we sample the zero moves from a Gaussian and we 
take no = 1. Standard Monte Carlo (l| is given by to = 1 
and ni = 1. Standard MTS |6] is recovered with m = 2. The 
early rejection method |2, §14.3.2] is recovered with to = 1 
(for all I). Note that a level-m move entails 111™=/ 1 n i' level-/ 
moves. At any stage in the recursion, we have the freedom to 
vary some of the components of E(T) adiabatically. 

In the following sections, we apply these techniques to 
constrained molecules. In simple cases, we can apply the 
MTS method semi-analytically to derive a correct constrained 
move. In more complicated cases, we apply the adiabatic 
technique to lift and to reapply the constraint. 

III. STIFF MOLECULES 

A constrained molecule is a mathematical idealization of a 
real system in which some degrees of freedom are stiff, i.e., 
the associated energies are large. Thus we can split the energy 
into "hard" (h) and "soft" (s) components, 

E{T) = E h (T) + E s (T), 

where T is the configuration of the system. For example, let us 
assume that an all-atom force field, such as Amber 01, pro- 
vides an accurate description of the system. (We recognize, 
of course, that present-day force fields are only approximate. 
However, our purpose here is to make the connection between 
an all-atom representation and a simpler rigid representation 
and, in this context, the details of the all-atom model are of 
secondary importance.) Then E^ might represent the bond 
stretching and bond bending terms, while E s is given by the 
other terms (bond torsion and the non-bonded energies). 

The constrained limit is now given by £j, — > 00. Before 
we consider this limit, it is useful to examine how the stiff 
system may be treated. Conventional Monte Carlo is ineffi- 
cient because, in order to have an reasonably large acceptance 
rate, the step-size needs to be set to a small value (determined 
by .Eh) so that diffusion in the soft directions is very slow. 
However, we can apply MTS Monte Carlo in this case with 
Ei = E h and E 2 = E s . 

Let us apply this method to a system of "rigid" molecules, 
e.g., water molecules, taking E^ to include the intra-molecular 
energies (responsible for maintaining the rigidity) and E s to 
include the inter-molecular energies. Suppose the level-0 
moves consist of symmetrically displacing the atoms in each 
molecule. The result of the ni level- 1 Monte Carlo steps 
will clearly be a symmetric, independent, and nearly rigid dis- 
placement (translation and orientation) of each molecule. This 
configuration is then subjected to Boltzmann acceptance with 
the inter-molecular energies. In this case, we can easily pass 
to the constrained limit (with exact rigidity), merely by ensur- 
ing that the trial (level- 1) moves of the molecules are rigid. 
In this case, we have just rederived the "standard" move for a 
system of rigid molecules. 

In order to illustrate the application to flexible molecules, 
we shall treat the molecules as being made up of several rigid 
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subunits or "fragments" connected by flexible bonds. How- 
ever we are interested in the limit where the inter-fragment 
bonds constrain the relative motions of fragments in certain 
ways, either by fixing the bond lengths (allowing the bond an- 
gles and bond dihedrals to vary) or by fixing the bond lengths 
and bond angles (allowing the bond dihedrals to vary). Such 
a model is adequate to describe a wide range of interesting or- 
ganic molecules including proteins and drug-like ligands. We 
assume that the rigidity of the fragments is imposed only by 
intra-fragment energy. If other terms (e.g., an improper tor- 
sion term involving atoms from two fragments) contribute to 
the rigidity of a fragment, then we shall treat such terms as 
additional inter-fragment energies. 

We apply the generalized MTS method to this system with 
m = 3, the intra-fragment energy given by E\, the inter- 
fragment bond constraints given by E2, and with E3 account- 
ing for all the other energies. The argument given above al- 
lows us to pass to the limit of strictly rigid fragments. The 
method is then equivalent to a standard MTS method where 
the "elementary" moves consist of rigid displacements of 
each fragment which are Boltzmann accepted with energy 
Eh = E%. A sequence of n = 712 such moves are made 
with the result Boltzmann accepted with energy E s = £3. A 
possible prescription |8l §VII] for the rigid displacements of 
the fragments is to translate the fragment by a vector sampled 
from an isotropic 3-dimensional Gaussian and to rotate the 
fragment by s about an axis s where s is a "rotation vector" 
also sampled from an isotropic 3-dimensional Gaussian. The 
variances for the two Gaussians should be adjusted so that the 
translational and rotational components result in comparable 
displacements of the atoms of the fragment. 

Provided that the inter-fragment constraint terms E^ are 
sufficiently stiff, it is not important to include a detailed 
model of these terms; because the motion will take place near 
the bottom of the constraint potential well, a harmonic (i.e., 
quadratic) approximation to the constraint potential will suf- 
fice. On the other hand, if the stiffness of the constraint energy 
depends on any of the soft degrees of freedom, it is important 
that this effect be included. 

It is frequently the case that Eh may be computed much 
more rapidly than E s . For example, when imposing bond 
constraints on a molecule, E^ requires O(N) computations, 
where N is the number of atoms, while E s requires 0(N 2 ) 
computations for the electrostatic and implicit solvation ener- 
gies. Thus we might be able to take n reasonably large and 
still have the computational cost dominated by the evaluation 
of E S (T). 

In order to realize the full benefits of imposing constraints 
we need to pass to the constrained limit (E^ — > 00). In 
this limit, the motion collapses onto a lower-dimensional sub- 
manifold in configuration space. Unfortunately, in contrast to 
the case of rigid molecules, we cannot appeal to symmetry to 
enable us to take this limit analytically. Instead, we use the 
adiabatic technique. 



IV. ADIABATICALLY VARYING THE STIFFNESS 

Let us rewrite the energy of the system, multiplying the 
E h (T) by T/T*, where T is the temperature of the system, 
and T* is a "constraint" temperature. The Boltzmann factor 
exp(— j3E), will then have the form 

exp(-/?£) = exp(-/3S s - P*E h ) 

where (3* = l/(kT*). 

In our application, where we are interested in the con- 
strained limit T* — > 0, a direct application of the MTS method 
leaves us with two bad choices. If we take T* to be sufficiently 
small that we can consider the constraints to be satisfied, we 
will have to chose the step size for the E^ moves to be so 
small that the change in configuration after n E^ moves will 
be small. On the other hand, letting T* be sufficiently large to 
allow moves will result in configurations where the constraints 
are poorly satisfied. 

We overcome this difficulty by regarding T* as a parameter 
(taking the place of A) and by adiabatically varying T* from 
zero (where the constraints are satisfied but MTS is ineffective 
at making moves) to a finite value (where the constraints are 
relaxed and MTS becomes effective) and back to zero again 
(to reimpose the constraints). During the course of changing 
T*, we make n E^ moves (each with the instantaneous value 
of T*). The effect of these n moves will be an E^ move with 
T* = 0, i.e., a move which satisfies the E^ constraint. 

It remains to give a recipe for varying T*. As we vary T*, 
we would naturally adjust the step size for the moves in such a 
way that the number of steps needed to equilibrate the system 
is a constant, suggesting that we vary T* exponentially. We 
therefore pick 



T * = [ TXexp(a(i - 1)), 
* I TXexp(a(n-i)), 



for < i < to, 
for m < i < n, 



where we have taken n = 2m + 1 and where T* is the con- 
straint temperature used for the zth E& move, T * = T* A is 
some temperature sufficiently small that we can consider the 
constraints to be exactly satisfied, and a is the rate of increase 
of the temperature which should be sufficiently small that the 
adiabatic condition is satisfied. Even though T\ and a are 
small, we can pick n sufficiently large that T^ +1 
T* A exp(am) is finite. 

In addition, we choose the step size for the ith move to 
be di = k^fT* where k is a constant. In traditional Monte 
Carlo, we normally pick k to maximize the diffusion rate 
which at the ith step is roughly 
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where A is the mean acceptance rate and (...) denotes an en- 
semble average. Maximizing the diffusion rate usually results 
in a rather small acceptance rate A ~ 0.1 because rare large 
steps can lead to faster diffusion than frequent small steps. 
However, in our application, where we want the system to 
remain in equilibrium as we vary the temperature, rare large 
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steps are bad. So we pick k to maximize ADi and this will 
usually result in A ~ 0.5. Note that for a given k, we have 

A ~ CT*, 

where C is constant provided that the step size is not too large. 
The overall diffusion can be estimated by summing over the n 
steps, 



D = 



<(r„-r ) 2 ) 



where we have assumed that successive steps are uncorrected 
and we have taken a <C 1 and Tg » Tj[. We should select 
parameters, a and Tjj, in order to adjust D so that the E s 
acceptance rate is 0(1). 

This method includes internal diagnostics to verify that a is 
small enough. We define . . .f (resp. . . .J.) as the average of a 
quantity over the steps where T* is increasing, i.e., i < m + 1 
(resp. decreasing, i.e., i > m + 1). We monitor A(ri)/T*f 
and £"h(ri)/T*| and demand that both should be close to the 
equilibrium value of N/2 (where N is the number of hard 
degrees of freedom). If a is too large, then we would find 



E h (TA/T*r « N/2, 



E h (Ti)/T t *i » N/2. 

In particular, if the final £^(r n ) is many times T* A , then the 
configuration is "hung up" and does not obey the constraints. 
If this happens frequently, the simulation needs to be rerun 
with a smaller setting for a; if, on the other hand, it happens 
only rarely, we would merely reject the step. We can also 
monitor the mean acceptance rates A~\ and A[. These should 
be about the same; however, if a is too large, we will find 
A] » A[. 

A useful guideline for picking TJ is that once the n E^ 
moves are completed and the system is presumably equili- 
brated to T\, we should be able to enforce the constraints by 
setting T* = (using any convenient energy minimization 
technique) with a negligible change in the configuration, e.g., 
with a negligible change in E S (T). 



V. PAIRWISE TERMS IN ENERGY 

Having made an adiabatic move using E^, the final step 
is to accept the move depending on the change in E s . We 
wish to compute this energy as efficiently as possible by us- 
ing the rigidity of the fragments. Force fields such as Amber 
01 include two types of energies: interactions between atoms 
(the electrostatic and Lennard- Jones terms) and bond energies 
(stretch, bend, and torsion). Since the number of terms in 
non-bonded energies typically scales as 0(N 2 ) where N is 
the total number of atoms in the system, while the number 
of bond terms scales as O(N), we concentrate on optimizing 
the evaluation of the non-bonded terms. In our case where 
the molecules consist of rigid fragments connected by flexi- 
ble bonds we need only include the bond terms contributed 



by the much smaller number of inter-fragment bonds. Fur- 
thermore, we need only include the energy contributed by the 
"free" components of such bonds. Thus, if the lengths and an- 
gles of such bonds are constrained, then we need only include 
the torsion energy in E S (T). 

We start by assuming that the non-bonded energy terms can 
be expressed as a sum over atom pairs. This applies to the 
electrostatic and Lennard-Jones terms in Amber H. However, 
implicit solvent models have a more complex structure and we 
consider these in the next section. 

Suppose our molecular system consists of N atoms. These 
atoms are grouped into M molecules and we denote Mi as the 
set of atoms making up the Zth molecule. Similarly, the atoms 
are divided into F rigid fragments and we denote F a as the 
set of atoms making up the ath fragment. A typical pairwise 
energy term can then be written as 

0<i<j<N 

where g denotes the type of energy term (electrostatic or 
Lennard-Jones), % and j are atom indices, nj is the distance 
between atoms i and j, f g is some function of distance, and 
C 9: ij is a coefficient which depends on the atoms but not on 
their positions. Thus for electrostatic interactions, Cgjj de- 
pends on the partial charges on the two atoms (assumed to be 
constant in Amber) and on the bonding relation between the 
atoms. Physical energy functions satisfy rim r _,oo f g (r) = 0. 
When the fragments are separated sufficiently, we have 



E„ 



E, 



go 



Cg,ijfg( r i 



Q<a<F i<i 



which is independent of T. It is convenient to choose E g g as 
the "origin" for the E g , i.e., we compute only 



E, 



gl 



E„ — E, 



go 



E 



E 

0<a<b<F i^Fa. 



We note that only energy differences enter into the computa- 
tion of observable quantities, and so we are free to select the 
arbitrary origin for energies. 

Let us consider the application of a small molecule (Ni 
atoms) interacting with a protein (N p ^> Ni atoms) where 
only some of the protein side chains near the binding site are 
allowed to move. By avoiding computing the interaction en- 
ergy between atoms in the immobile portion of the protein, 
the above prescription reduces the computational cost from 
0(N*) to 0(N t N p ). 

This cost may still be too large and we can substantially 
reduce the cost by implementing energy cutoffs for the inter- 
actions. This is easily accomplished by multiplying f g (rij) 
by a cutoff function, c g (rij). A possible form for this cutoff 
function is 



i( r gi) 



r g2 



r g2 - r g i 



forr < r gl , 
forr > r g2 , 

otherwise, 
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with r g i < r g 2, which linearly tapers the energy to zero over 
[r g i, r g2 ). Other tapering functions can be employed, or, by 
choosing r g 2 = r g \, we can implement a sharp cutoff. This 
type of cutoff function implements a per-atom cutoff and is 
appropriate for energy terms which are additive at large dis- 
tances, such as the Lennard- Jones potential. The electrostatic 
potential, however, involves substantial cancellation at large 
distances — two neutral molecules interact via a dipole-dipole 
term which varies as 1/r 3 , while the individual atom-atom 
terms decay as 1/r.y. In this case, we need to identify groups 
of atoms which should interact together. The residues of a 
protein provide a convenient grouping and we would typically 
assign all the atoms in a small-molecule ligand to a single 
group. Compatible with the usage for a protein, we refer to 
these groups as residues. For each residue, s, we define a cen- 
ter position, h s , most conveniently defined as the center of 
mass, and a radius, defined as the radius h s of the sphere cen- 
tered at b s which includes the van-der-Waals spheres of radius 
Pi of all the constituent atoms. We then apply a "per-residue" 
cutoff function multiplying the contribution from the residue 
pair (s, t) by c 9 (|b s - b t \-(h s + h t )). 

The values used for the cutoff radii, r g \ and r g 2, need to 
evaluated based on the accuracy desired for the simulation. 
This can be determined by numerically determining the dif- 
ference in the results (either for the energies directly or for 
some derived quantity such as binding affinity) between the 
finite- and infinite-cutoff energies. In applications to Monte 
Carlo codes, it is possible to carry out the sampling at an en- 
ergy approximating the actual energy and to compensate for 
this when performing the canonical averages (which might be 
carried out on a subset of the Markov chain). In this case, 
the sampling energy might entail using shorter cutoffs than 
would be warranted on the basis of accuracy. Having deter- 
mined suitable cutoffs, it is a simple matter to evaluate the 
energy avoiding treating atom pairs beyond the respective cut- 
offs. In the following, we treat electrostatic (e) interactions, 
with a per-residue cutoff, and Lennard- Jones (I) interactions, 
with a per-atom cutoff; furthermore we assume that r e2 > ri 2 , 
i.e., the electrostatic interactions are longer range than the 
Lennard- Jones. 

We first loop over all the atoms in each residue computing 
b s and h s for all residues s. We then loop over all pairs of 
residues, s < t, skipping any pair whose atoms all belong 
to the same fragment or those for which |b s — b f | > r e2 + 
h s + h t . If the residue pair survives these tests, then all atom 
pairs from different fragments are considered; if s = t, 
we restrict the pairs to i < j. All such pairs contribute to 
the electrostatic energy while those which satisfy < 
contribute to the Lennard-Jones energy. There obviously is 
scope for additional optimization here. For example, the inner 
atom loop can be skipped if the second residue belongs to a 
single fragment which matches the fragment of a particular 
atom in the first residue. 

Because of the way in which the cutoffs are applied, the re- 
sult for the energy is independent of the assignment of atoms 
to residues for energy terms which use a per-atom cutoff. In 
addition, differences in the non-bonded energies are indepen- 
dent of the assignment of atoms to fragments. The energies 



for assemblies of 3 or more molecules can be expressed in 
terms of the energies of 1 or 2 molecules. These provide use- 
ful checks on the implementation. 

In some contexts it is useful also to define a "steric" energy 
term which is infinite if any atoms overlap (with some defini- 
tion of a "hard" atom radius) and is zero otherwise. This pro- 
vides a rapid check of new configurations — particularly when 
trying to "insert" a molecule during a grand canonical sim- 
ulation 1911 or when switching systems using the wormhole 
method lHoll . A conservative definition of the hard atom ra- 
dius is 0.55pi for non-bonded atom pairs and 0A5pi for 1-4 
atom pairs. We skip the check for 1-2 and 1-3 pairs and for 
those atoms with pi — 0. This energy term can be imple- 
mented in essentially the same way as described above but 
with scope for additional speedups. The cutoff radius in the 
residue-residue distance check can be replaced by 0. An ad- 
ditional atom-residue distance check can be be used to avoid 
executing the inner atom loop if the outer atom is outside the 
sphere for the second residue. Finally, as soon as an overlap 
of hard spheres is detected the routine can immediately return 
an infinite result. 



VI. IMPLICIT SOLVENT MODELS 

We now turn to the computation of the energy term for im- 
plicit solvent models. We focus here on the generalized Born 
solvent models full an d we have considered various imple- 
mentations QEHEEllI- Evaluating the solvation 
energy for a system of molecules with such models is typi- 
cally orders of magnitude slower than computing the energy 
of the molecules in vacuum. The computation time is fre- 
quently compared to the time to compute the energy with an 
explicit solvent model (including O(10 3 ) solvent molecules). 
However, such comparisons are misleading because implicit 
solvent models do not attempt to compute the energy of a 
particular configuration of solvent molecules but to compute 
the free energy of solvation, i.e., to average over all possi- 
ble solvent configurations for a given configurations of solute 
molecules. Thus the chief benefit of an implicit solvent model 
is to reduce dramatically the number of degrees of freedom 
in the problem. In the generalized Born solvent models, the 
energy is written as the sum of two terms: a polar term which 
is usually called the "GB" term and a cavity term which is 
proportional to the solvent accessible surface area, the "SA" 
term. 

The GB term involves long-range interactions and is the 
most costly to compute. We address the calculation of this 
term first. The basic expression is 111 ill 

Gpoi = ~\^r\ l ~ ~) ^2<liQjf{rij,ai,aj), (1) 

where e s is the permittivity of the solvent, f(nj,ai,aj) = 
[rfj + OLiOt.^ ■ exp(— r| J -/(4a i aj))] -1 / 2 , and the double sum 
runs over all pairs of atoms (including i = j and i «s j). 
In eq. {!}, a, is the "generalized" Born radius of the ith atom, 
which is larger that the "bare" Born radius to account for the 
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fact that atoms close to i partially shield it from the solvent. 
Gp i represents the electrostatic energy required to solvate a 
pre-assembled group of molecules and thus this term is added 
to the vacuum electrostatic energy. The various implementa- 
tions for the GB term differ in how oti is computed. 

For illustrati ve purp oses, let us consider the model of 
Hawkins et al. (OF 13 E (With minor modifications, 
the technique is applicable to other GB models.) We express 
Oi as eq. (10)] 



1 1 A 



where pi is the radius of atom i, 

dr 



A, 



Hij (r; Tij , pj) 



(2) 



(3) 



is the reduction in the effective inverse Born radius of atom i 
due to atom j. Here Hij is the fraction of the area of a sphere 
of radius r centered on the ith atom eclipsed by a jth atom 
and is given by fl3l eq. (12)] 



Pi - in 



H 



1. 
0, 



ArijT 



for |r 



Pj ' '' ' r ij +Pj, 
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zero for pj > pi + r^, which is a possibility not considered 



in|13J. 

Clearly G po \ is no longer the sum of pairwise atom-atom 
contributions because the interaction of two atoms is affected 
by the modification of the dielectric environment by a third 
atom. However G po i may be evaluated by two pair-wise op- 
erations carried out in sequence. The first evaluates the gen- 
eralized Born radii a, and the second computes the resulting 
electrostatic energy. 

As with the treatment of the electrostatic and Lennard- 
Jones terms, we can seek to limit the computational cost 
of evaluating G po i by the use of cutoff functions. Because 
eq. Q provides the dielectric screening for the vacuum elec- 
trostatic term, it is important that the cutoff function multiply- 
ing f(rij, ai, aj) exactly match that used for the electrostatic 
term. 

We also introduce a cutoff in eq. (0 by multiplying Ay by 
Cb(rij). A per-atom cutoff is justified since all the Ajj are 
positive. Because Aij scales as 7y for large r v -, the error 
introduced by Cb(rij ) scales relatively slowly as rZ, . In prac- 
tice, this means we need to make r^i reasonably large which 



in turn means that the cost of evaluating G po i in the case of a 
small ligand interacting with a protein is much larger than the 
cost for the electrostatic potential. In particular, the screening 
of the ligand may modify the Born radii of a large number of 
protein atoms and this unavoidably leads to a large number of 
pair contributions to eq. 0. 

The procedure for computing the energy outlined in the pre- 
vious section can now be modified to deal with the evaluation 
of Gp i. As before our "zero" energy is given by separat- 
ing all the fragments of all the molecules infinitely far apart. 
We set up the calculation of a system of molecules by pre- 
computing a-to which is given by eq. (0 with the sum restrict- 
ing to include only the intra-fragment contributions (i.e., index 
j ranges only over atoms within the same fragment as atom i). 
We compute Ay and Aji together because they involve many 
of the same terms, allowing the loops to be restricted to i < j, 
and we apply the Born cutoff to the calculation of o^o . 

When computing the energy of a molecular system, we 
compute all the updates to the Born radii due to atoms in dif- 
ferent fragments within the Born cutoff, applying the same 
techniques of lumping the atoms into residues described above 
(which allows the cutoff criteria to be applied to groups of 
atoms) and of restricting the loops to s < t and, for s = t, 
to i < j. During this phase we mark all the residues which 
contain atoms with ^ a;o- We then make a second pass 
over the atoms to evaluate the terms in eq. Q. We use the 
i ^ j symmetry of the summand to make the restrictions 
s < t and, for s = t, i < j. In the innermost loop, we accumu- 
late qiqjf(rij , ai, aj) if i and j belong to different fragments. 
Otherwise, we add (/,</, ./'•:/•,,.'.,. aj) - f(rij,a i0 , a j0 )] and 
we can skip this evaluation if both ai — a^ and aj = ajQ. In 
addition, we can skip pairs of residues if all the atoms in each 
residue belong to the same fragment and if neither residue is 
marked as having modified Born radii. 

Salt effects 1 18] are easy to include within this framework. 
A minor complication occurs in the GB model of Qiu et 
al. 81211 because ct;o depends on the "volume" of the atoms and 
in this model the volume depends on the 1-2 bonded atoms 
which may belong to a different fragment. We account for 
this by assuming the presence of such bonded atoms with an 
ideal bond length. This is, therefore, only exact if the inter- 
fragment bonds are at their ideal lengths. Our treatment here 
may be considered as a generalization of the frozen atom ap- 
proximation for GB/SA 1 19]. However, in our application we 
make all the approximations in the energy function and the re- 
sulting energy is then a "state variable" and simulations based 
on this are well behaved. In contrast the implementation of 
frozen atom approximation defines the energy so that it de- 
pends on the history of the system which may cause the sim- 
ulation to exhibit unphysical properties. 



VII. SOLVENT ACCESSIBLE SURFACE AREA 

The other important contribution to the solvation free en- 
ergy is the cavity term. This is obtained by placing spheres 
centered at each atom with radius a\. = pi + r w where r w is 
a nominal water radius (typically r w = 0.14 nm). The cavity 
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term is given by 

^cav ^ ^ ® jA-j , 

where is the "solvent accessible surface area" for the ith 
atom, i.e., the exposed surface area of the spheres around i 
which is not occluded by any other spheres and oi is the sur- 
face tension for the ith atom. (Typically oi is taken to be a 
constant independent of atom, Oi ss 3 kJ mol -1 nm~ 2 ; how- 
ever the method we describe does not require this assump- 
tion.) As before, the zero energy state is obtained by sep- 
arating the fragments infinitely. The energy is then given 
by the additional occlusion of the surface that occurs as the 
fragments are assembled into molecules and the molecules 
brought into contact with one another. 

The exact evaluation of this term is quite complex and for 
this reason a simple pairwise approximation has been devel- 
oped 1 20]. However, the errors in this method are poorly 
quantified. This together with the fact that this term is typi- 
cally small compared to the electrostatic terms in the energy 
lead us to develop a simple zeroth-order quadrature method. 
We select an accuracy level for the cavity calculation 8, e.g., 
S = O.lkJ/mol. We prepare for the calculation of the 
cavity term by placing each fragment in a "template" posi- 
tion and we arrange a set of points on a sphere of radius a; 
around each atom i. The number of points is chosen to be 
Ni = \\-Ka\uijS\. The points are distributed approximately 
uniformly around each sphere and the entire surface energy 
of the sphere, Airajdi is divided among the Ni points. (We 
will discuss the details of how to select the points and assign 
the energy later.) We next perform the intra-fragment occlu- 
sion by deleting all the points of atom i which are within a,j of 
some atom j ^ i. In this way each fragment is surrounded by 
a cloud of surface points each representing about 6 of cavity 
energy. 

In order to compute the cavity term for a particular molec- 
ular configuration we transform the surface points for each 
fragment from their template positions to their actual positions 
and make a copy of the cavity energies for each point. We con- 
sider all pairs of atoms (i, j) such that i and j are in different 
fragments and nj < a j + aj . We subtract from G cav the ener- 
gies of all the points on atom i that are within aj of atom i and 
we set the energies of these points to zero (to avoid their being 
counted multiple times). The optimizations described above 
can be used: the application of a residue-residue cutoff (ex- 
cluding residue pairs (s,t) with |b s — b t | > 2r w +h s +h t ), an 
atom-residue cutoff, and the treatment of the and 
terms together. 

In practice, the cost of evaluating this term is small for 
5 w O.lkJ/mol. The error is proportional to S and it is 
easy to benchmark a particular calculation by repeating it with 
smaller S. The resulting G cav is obviously a discontinuous 
function of configuration, jumping by ±<5 as points move in 
and out of the water spheres of other atoms. Thus it's an inap- 
propriate model for a molecular dynamics simulation. How- 
ever, it yields satisfactory results for Monte Carlo simulations. 

Let us return to the question of how to position the points 
on the atom sphere and how to divide the energy between 



these points. Ideally, we would divide the energy of the sphere 
based on the area of Voronoi polygons around each point. The 
error will then be proportional to the maximum radius of the 
Voronoi polygons and the ideal distribution of points is the one 
which minimizes this maximum radius. This is the so-called 
"covering problem" for the sphere, i.e., how to cover a sphere 
with identical discs ll2lll . Unfortunately, there are no general 
solutions to this problem. So instead we divide the sphere into 
equal intervals of latitude and we divide each latitudinal inter- 
val longitudinally into approximately square regions. A point 
is placed at the center of each region and the area of the re- 
gion is assigned to that point. Within each fragment, we alter 
the position of the pole from one atom to the next, in order 
to avoid the occlusion of many points simultaneously as frag- 
ments move relative to one another. 



VIII. DISCUSSION 

We have shown how to make Monte Carlo moves for a 
molecular system with constraints. Constraints are imposed 
in a realistic way ensuring that we obtain the right distribu- 
tion corresponding to a thermodynamic equilibrium. We will 
still need to know this constrained distribution if we wish to 
make wormhole moves fioll . because, in order to satisfy de- 
tailed balance, we require knowledge of the wormhole vol- 
umes and these include a factor proportional to the "thick- 
ness" of the constraint manifold. The adiabatic move involves, 
naturally, many evaluations of the constraint energy raising a 
concern that the implementation will be slow. In reality, the 
cost of evaluating the constraint energy is minuscule, particu- 
larly in comparison with the solvation energy, so it is possible 
to evaluate the constraint energy many thousands of times in 
the course of an adiabatic move with minimal impact on the 
overall running time. The method avoids much of the algebra 
associated with other ways of imposing constraints 1 221 and 
thus is more flexible and is easier to implement. 

In the simple case of a molecule in which only a number 
of dihedral angles are allowed to vary, the movement of all 
the atoms in the molecule is bounded and thus the soft-energy 
acceptance probability is reasonably large. In contrast, the 
method where the dihedral angles are perturbed may lead, due 
to a lever effect, to large motions if the molecule itself is large. 

This method can easily be generalized to do localized 
movements. Thus, we can tailor the random displacements 
of a protein to explore the movement of a single loop. De- 
tailed balance is ensured if the random displacement is a func- 
tion of the atom but not of its position. (The general case 
can be accommodated by a suitable factor in the acceptance 
probability.) This method of localized movements is more 
widely applicable than techniques such as "concerted rota- 
tions" I22i l23l 12411 . Artificially fixing the positions of some 
atoms would, of course, mean that the moves would not be 
ergodic. This would be justified if we were interested in ex- 
amining the restricted system and we would then require er- 
godicity over the restricted configuration space. 

We have also considered how to optimize the evaluation 
of the energy in a system of molecules made up of rigid 
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fragments bonded together. This allows the use of implicit 
solvent at an acceptable cost. If the system is further con- 
strained to allow only the variation of the torsion angle of the 
inter-fragment bonds (fixing the bond lengths and bond an- 
gles), then we should also consider modifying the force field 
to "loosen" the torsion energies to counteract the effect of 
the hard constraints on the other bond terms. G5 and Scher- 
aga 1 25] show the importance of considering such an effect 
and Katrich et al. ll2al have offered a prescription for con- 
verting a general force field to include this effect. Alterna- 
tively, we might consider re-parameterizing the torsion terms 
by carrying out constrained geometry optimizations of model 
molecules where the energy of the molecule is minimized with 
the dihedral angles fixed 1 27] . 
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